## ----setup, include=FALSE------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, message=FALSE, warning=FALSE)
library(Zelig)
library(MASS)
library(BayesTree)
library(birdring)
library(dplyr)
library(lmtest)
library(car)
library(ggplot2)
effects <- coefs <-  vector(mode="list", length=3)

## ------------------------------------------------------------------------
setwd("~/Dropbox/Zach/JLC_replication")
load("cert_data.rda")
meanSD <- function(x, digits=3){
   out <- c(mean(x, na.rm=T), sd(x, na.rm=T))
   out <- round(out, digits)
   names(out) <- c("Mean", "SD")
   out
}
shuffle <- function(x,y){
   res <- matrix(c(rbind(x,y)), ncol=1)
   colnames(res) <- "B (HPD)"
   rownames(res)[seq(2, nrow(res), by=2)] <- ""
   rownames(res)[seq(1, nrow(res), by=2)] <- names(x)
   noquote(res)
}
mci <- function(x)c(mean(x, na.rm=T), quantile(x, c(.025, .975), na.rm=T))
x2test <- function(L, rhs, b, V){
    value.hyp <- L %*% b - rhs
    vcov.hyp <- L %*% V %*% t(L)
    SSH <- as.vector(t(value.hyp) %*% solve(vcov.hyp) %*% value.hyp)
    q <- NROW(L)
    p <- pchisq(SSH, q, lower.tail = FALSE)
    return(p)
}

## ------------------------------------------------------------------------
x <- as.matrix(cert_data[,c("casecert", "lctdir", "lowercourtdissent", "abortion", "deathpenalty", "allegationconflict", "constitutionalclaim", "lowercourtreversal", "uspetitioner", "trueconflict", "amicusfavor", "amicusagainst")])
x <- sapply(1:ncol(x), function(z)as.numeric(x[,z]))
cc <- complete.cases(x)
y <- x[which(cc),1]
x <- x[which(cc),-1]
tab <- round(t(do.call(rbind, by(x, list(y), colMeans))), 2)
ntab <- c("Lower Court Direction", "Lower Court Dissent", "Abortion", "Death Penalty", "Alleged Conflict", "Constitutional Claim", 
   "Lower Court Reversal", "US Petitioner", "True Conflict", "Amicus Favor Cert", "Amicus Against Cert")
rownames(tab) <- ntab
round(tab, 2)

## ------------------------------------------------------------------------
x <- as.matrix(cert_data[,c("casecert", "lctdir", "lowercourtdissent", "abortion", "deathpenalty", "allegationconflict", "constitutionalclaim", "lowercourtreversal", "uspetitioner", "trueconflict", "amicusfavor", "amicusagainst")])
x <- sapply(1:ncol(x), function(z)as.numeric(x[,z]))
cc <- complete.cases(x)
y <- x[which(cc),1]
x <- x[which(cc),-1]

## ------------------------------------------------------------------------
# estimate model 
pmod <- glm(y ~ x, family=binomial)
# generate the posterior draws of coefficients from the logistic regression model 
set.seed(1256)
b.parm <- mvrnorm(1500, coef(pmod), vcov(pmod))
# generate posterior draws of cert-worthiness
parm.post <- plogis(model.matrix(pmod) %*% t(b.parm))
# transpose matrix of posteriors to match output from BART
parm.post <- t(parm.post)

## ------------------------------------------------------------------------
# Estimate BART
set.seed(1120477)
out <- bart(x, y, ndpost=1500)
# Calculate posterior probabilities from fitted values
bart.post <- pnorm(out$yhat.train)

## ---- fig.show='hide'----------------------------------------------------
# initialize the object to hold the overlap coefficients
olaps <- rep(NA, ncol(bart.post))
# calculate the overlap for each column in the matrices of posteriors
for(i in 1:length(olaps)){
    olaps[i] <- overlap(parm.post[,i], bart.post[,i], edge.of.parameter.space=T)
}
# identify the outlying cases
o <- which(colMeans(parm.post) < .55 & colMeans(bart.post) > .65)
# summarize all overlaps
summary(olaps)
# summarize overlaps for the non-outlying cases
summary(olaps[-o])
# summarize overlaps for only outlying cases. 
summary(olaps[o])

## ------------------------------------------------------------------------
# calculate difference between two matrices of probabilities
diffs <- parm.post - bart.post
# for each column, figure out what proportion of those differences are bigger than 0
pv <- apply(diffs, 2, function(x)mean(x > 0))
## Note no p-values outside the range (0.05, 0.95)
summary(pv)

## ---- fig.height=8, fig.width=8, out.width=".65\\textwidth", out.height=".65\\textwidth", fig.align="center"----
## Figure A1
tiff("figureA1.tiff", width=8, height=8, units="in", res=300)
plot(colMeans(parm.post), colMeans(bart.post), xlab="Logistic Regression Posterior Mean Probabilities", ylab="BART Posterior Mean Probabilities")
ol <- which(colMeans(parm.post) < .55 & colMeans(bart.post) > .65)
ol0 <- ol[-8]
ol1 <- ol[8]
text(colMeans(parm.post)[ol0], colMeans(bart.post)[ol0], cert_data$docket[ol0], 
    pos=c(4, 2, 2, 2, 2, 2, 2,  2, 2, 2, 4))
text(.556, .903, cert_data$docket[ol1], 
    pos=1)
dev.off()

## ------------------------------------------------------------------------
cert_probs <- array(dim=c(nrow(cert_data), 1500))
cert_probs[which(cc), ] <- t(parm.post)
cert_probs <- as.data.frame(cert_probs)
cert_probs$docket <- cert_data$docket

## ------------------------------------------------------------------------
cert_post <-  left_join(data.frame(docket=moddat$docket), cert_probs)
cert_post <- cert_post[,-1]
moddat$certworthy1 <- rowMeans(cert_post)

## ------------------------------------------------------------------------
bv <- moddat$blackmunvoteoncert2
mc <- moddat$mcrecommendation
bc <- moddat$blackrecommendation
bv <- ifelse(is.na(bv), 0, bv)
mc <- ifelse(is.na(mc), 0, mc)
bc <- ifelse(is.na(bc), 0, bc)
temp <- data.frame(bv=bv, bc=bc, mc=mc, casecert=moddat$casecert)

## ------------------------------------------------------------------------
prop.pop <- c(.85, .15)
t.mc.mod <- glm(mc ~ as.factor(casecert)-1, data=temp, family=binomial)
set.seed(64502)
p.mc <- plogis(mvrnorm(1500, t.mc.mod$coef, vcov(t.mc.mod)))
t.mc <- p.mc %*% prop.pop
meanSD(t.mc)

## ------------------------------------------------------------------------
# Blackmun's Clerk
t.bc.mod <- glm(bc ~ as.factor(casecert)-1, data=temp, family=binomial)
set.seed(1205604)
p.bc <- plogis(mvrnorm(1500, t.bc.mod$coef, vcov(t.bc.mod)))
t.bc <- p.bc %*% prop.pop
meanSD(t.bc)

# Blackmun
t.bv.mod <- glm(bv ~ as.factor(casecert)-1, data=temp, family=binomial)
set.seed(1105042)
p.bv <- plogis(mvrnorm(1500, t.bv.mod$coef, vcov(t.bv.mod)))
t.bv <- p.bv %*% prop.pop
meanSD(t.bv)

## ------------------------------------------------------------------------
tau <- data.frame(mc = t.mc, bc=t.bc, bv=t.bv)

## ------------------------------------------------------------------------
black.ideo <- c( -.7661515 , -.9495762 , -.9878716 , -1.001301 , -.8643821 , -1.15319 , -1.383177 , -1.545243 , -1.802584 )

## ------------------------------------------------------------------------
year <- as.numeric(gsub("(^\\d{2}).*", "\\1", moddat$docket))-85
moddat$black.ideo <- black.ideo[year]

## ------------------------------------------------------------------------
moddat$dist_black_mc <- abs(moddat$black.ideo-moddat$mqtermprior)

## ------------------------------------------------------------------------
w <- which(complete.cases(moddat) & complete.cases(cert_post))
cert_post <- cert_post[w,] 
moddat <- moddat[w,] 

b1 <- array(dim=c(1500, 4))
set.seed(12294292)
for(i in 1:1500){
    m1a <- relogit(mcrecommendation ~ lctdir + cert_post[,i] + mcdist, data=moddat, tau=tau$mc[i])
    b1[i, ] <- mvrnorm(1, coef(m1a), vcov(m1a))
}

## ------------------------------------------------------------------------
h1 <- apply(b1, 2, coda:::HPDinterval.mcmc)
coef1 <- apply(b1, 2, mean)
p1 <- apply(b1, 2, function(x)mean(x > 0))
ph1 <- paste("(", apply(round(h1, 2), 2, paste, collapse=", "), ")", sep="")
pb1 <- round(coef1, 2)
names(pb1) <- c("Intercept", "Lower Court Direction", "Cert Worthiness", "Distance to Court Median")
shuffle(pb1, ph1)
coefs[[1]] <- shuffle(pb1, ph1)

## ------------------------------------------------------------------------
b1m <- colMeans(b1)
m1a <- relogit(mcrecommendation ~ lctdir + plogis(certworthy1) + mcdist, data=moddat, tau=mean(tau$mc))

## ------------------------------------------------------------------------
p1m <- plogis(model.matrix(m1a) %*% b1m)

## ------------------------------------------------------------------------
res <- NULL
gs <- seq(0.01, 0.99, length=1000)
for(x in gs){
    res <- c(res, mean((p1m > x) ==  model.response(model.frame(m1a))))
}

## ------------------------------------------------------------------------
y1 <- model.response(model.frame(m1a))
my1 <- mean(y1)
pm1 <- max(my1, 1-my1)
pre1 <- (max(res) - pm1)/(1-pm1)

## ------------------------------------------------------------------------
m1a <- relogit(mcrecommendation ~ lctdir + certworthy1 + mcdist, data=moddat, tau=tau$mc[i])
X1 <- X2 <- model.matrix(m1a)
X1[,2] <- 0
X2[,2] <- 1

## ------------------------------------------------------------------------
mes <- plogis(X2 %*% t(b1)) - plogis(X1 %*% t(b1))
bds <- coda:::HPDinterval.mcmc(colMeans(mes))
bds <- paste("(", round(bds[1], 3), ", ", round(bds[2], 3), ")", sep="")

## ------------------------------------------------------------------------
X1 <- X2 <- model.matrix(m1a)
X1[,3] <- 0
X2[,3] <- 1
mes2 <- plogis(X2 %*% t(b1)) - plogis(X1 %*% t(b1))
bds2 <- coda:::HPDinterval.mcmc(colMeans(mes2))
bds2 <- paste("(", round(bds2[1], 3), ", ", round(bds2[2], 3), ")", sep="")

## ------------------------------------------------------------------------
X1 <- X2 <- model.matrix(m1a)
X1[,4] <- 0
X2[,4] <- 5
mes3 <- plogis(X2 %*% t(b1)) - plogis(X1 %*% t(b1))
bds3 <- coda:::HPDinterval.mcmc(colMeans(mes3))
bds3 <- paste("(", round(bds3[1], 3), ", ", round(bds3[2], 3), ")", sep="")

## ------------------------------------------------------------------------
b2 <- array(dim=c(1500, 7))
set.seed(692042)
    for(i in 1:1500){
    m2a <- relogit(blackrecommendation ~ blackmundistance + mcrecommendation*dist_black_mc + cert_post[,i] + lctdir, data=moddat, tau=tau$bc[i])
    b2[i, ] <- mvrnorm(1, coef(m2a), vcov(m2a))
}

## ------------------------------------------------------------------------
h2 <- apply(b2, 2, coda:::HPDinterval.mcmc)
coef2 <- apply(b2, 2, mean)
p1 <- apply(b2, 2, function(x)mean(x > 0))
ph2 <- paste("(", apply(round(h2, 2), 2, paste, collapse=", "), ")", sep="")
pb2 <- round(coef2, 2)
names(pb2) <- c("Intercept", "Blackmun Distance to Median", "Memo Clerk Recommendation", "Blackmun Distance to Memo Clerk", "Cert-worthiness", "Lower Court Direction", "Memo Clerk Recommendation x Blackmun Distance to Memo Clerk")
shuffle(pb2, ph2)
coefs[[2]] <- shuffle(pb2, ph2)

## ------------------------------------------------------------------------
b2m <- colMeans(b2)
m2a <- relogit(blackrecommendation ~ blackmundistance + mcrecommendation*dist_black_mc + plogis(certworthy1) + lctdir, data=moddat, tau=mean(c(.072,.132)))
p2m <- plogis(model.matrix(m2a) %*% b2m)
res <- NULL
for(x in seq(0.01, 0.99, length=1000)){
    res <- c(res, mean((p2m > x) ==  model.response(model.frame(m2a))))
}
y2 <- model.response(model.frame(m2a))
my2 <- mean(y2)
pm2 <- max(my2, 1-my2)
pre2 <- (max(res) - pm2)/(1-pm2)

## ------------------------------------------------------------------------
m2a <- relogit(blackrecommendation ~ blackmundistance + mcrecommendation*dist_black_mc + certworthy1 + lctdir, data=moddat, tau=mean(tau$bc))
X2ma <- X2mb <- model.matrix(m2a)
X2ma[,6] <- 0
X2mb[,6] <- 1

mes <- plogis(X2mb %*% t(b2)) - plogis(X2ma %*% t(b2))
bds <- coda:::HPDinterval.mcmc(colMeans(mes))
bds <- paste("(", round(bds[1], 3), ", ", round(bds[2], 3), ")", sep="")

## ---- fig.height=8, fig.width=8, out.width="\\textwidth", out.height="\\textwidth", fig.align="center"----
# Set evaluation points for cert-worthiness
c2 <- seq(0, 1, length=4)
# generate hypothetical data for memo-clerk recommendation, distance of blackmun to memo clerk and cert-worthiness
eg <- expand.grid(mc = c(0,1), dist_black_mc =seq(.05, 4.55, length=10), certworthy1=c2)
# loop over the rows in the hypothetical data
res1 <- NULL
for(i in 1:nrow(eg)){
    #  make a copy of the data
    # set the variables of interest (mc recommendation, cert-worthiness and distance of Blackmun to mc
    # in our copy of the data to the values in the fake data for each row in turn
    tmp <- moddat
    tmp$mcrecommendation <- eg[i,1]
    tmp$dist_black_mc <- eg[i,2]
    tmp$certworthy1 <- eg[i,3]
    # make the design matrix for the new hypothetical data
    X1 <- model.matrix(formula(m2a), data=tmp)
    # replace cert-worthiness with the correct value
    X1[,5] <- eg[i,3]
    # calculate the posterior distribution of the predicted probabilities
    pp1 <- plogis(X1 %*% t(b2))
    # find the mean and confidence intervals for the predicted probabilities
    q1 <- mci(colMeans(pp1))
    # collect results
    res1 <- rbind(res1, q1)
}
# make results into a data frame with the hypothetical data, too.
colnames(res1) <- c("mean", "lower", "upper")
rownames(eg) <- rownames(res1) <- NULL
plot.dat <- cbind(eg, res1)
# change the values of cert-worthy in our plot data to a factor
plot.dat$certworthy1 <- factor(plot.dat$certworthy1, levels=sort(unique(plot.dat$certworthy1)), labels=paste("Pr(Cert=", round(sort(unique(plot.dat$certworthy1)), 3), ")", sep=""))
# change memo clerk's recommendation into a factor
plot.dat$mc <- factor(plot.dat$mc, labels=c("No", "Yes"))
# make the figure
# tiff("figure1.tiff", width=8, height=8, units="in", res=300)
mp <- ggplot(plot.dat, aes(x=dist_black_mc, y=mean)) + theme_bw()
mp + geom_ribbon(aes(ymin=lower, ymax=upper, fill=mc), alpha=.25) + scale_fill_manual("Memo Clerk\nRecommenation",values=c("gray25", "gray75")) + facet_wrap(~certworthy1) +
    geom_line(aes(x=dist_black_mc, y=mean, color=mc), show.legend=F) + scale_colour_manual(values=c("black", "black")) +
    guides(fill=guide_legend(title="Memo Clerk\nRecommendation"), title.position="top") + theme(legend.position="top") +
    labs(x = "Ideological Distance to Blackmun", y="Pr(Blackmun's Clerk Recommends Cert)")
# dev.off()

## ------------------------------------------------------------------------
X2ma <- X2mb <- model.matrix(m2a)
X2ma[,5] <- 0
X2mb[,5] <- 1

mes <- plogis(X2mb %*% t(b2)) - plogis(X2ma %*% t(b2))
bds <- coda:::HPDinterval.mcmc(colMeans(mes))
bds <- paste("(", round(bds[1], 3), ", ", round(bds[2], 3), ")", sep="")

## ------------------------------------------------------------------------
m3a <- relogit(blackmunvoteoncert2 ~ blackmundistance +  blackrecommendation*mcrecommendation*dist_black_mc + plogis(certworthy1) + lctdir,
    data=moddat, tau=mean(tau$bv))
m3b <- relogit(blackmunvoteoncert2 ~ blackmundistance + blackrecommendation +mcrecommendation*dist_black_mc +  plogis(certworthy1) + lctdir, data=moddat, tau=mean(tau$bv))
lrtest(m3a, m3b)

## ------------------------------------------------------------------------
n <- colnames(model.matrix(m3a))
linearHypothesis(m3a, paste(n[c(8:11)], "=0", sep=""))

## ------------------------------------------------------------------------
b3 <- array(dim=c(1500, 11))
set.seed(1561)
for(i in 1:1500){
    m3a <- relogit(blackmunvoteoncert2 ~ blackmundistance + blackrecommendation*mcrecommendation*dist_black_mc +
        cert_post[,i] + lctdir, data=moddat, tau=tau$bv[i])
    b3[i, ] <- mvrnorm(1, coef(m3a), vcov(m3a))
}

## ------------------------------------------------------------------------
L <- matrix(0, ncol=ncol(b3), nrow=4)
L[cbind(1:4, 8:11)] <- 1
rhs <- rep(0, 4)
x2test(L, rhs, colMeans(b3), var(b3))

## ------------------------------------------------------------------------
b3a <- array(dim=c(1500, ncol(model.matrix(m3a))))
b3b <- array(dim=c(1500, ncol(model.matrix(m3b))))
p3a <- p3b <- array(dim=c(1500, nrow(moddat)))
deva <- devb <- rep(NA, 1500)
set.seed(6603)
for(i in 1:1500){
    m3a <- relogit(blackmunvoteoncert2 ~ blackmundistance + blackrecommendation*mcrecommendation*dist_black_mc +
        cert_post[,i] + lctdir, data=moddat, tau=tau$bv[i])
    m3b <- relogit(blackmunvoteoncert2 ~ blackmundistance + blackrecommendation + mcrecommendation*dist_black_mc +
        cert_post[,i] + lctdir, data=moddat, tau=tau$bv[i])
b3a[i,] <- mvrnorm(1, coef(m3a), vcov(m3a))
b3b[i,] <- mvrnorm(1, coef(m3b), vcov(m3b))
p3a[i,] <- predict(m3a, type="response")
p3b[i,] <- predict(m3b, type="response")
deva[i] <- deviance(m3a)
devb[i] <- deviance(m3b)
}
y <- model.response(model.frame(m3a))
illa <- apply(p3a, 1, function(x)ifelse(y==1, log(x), log(1-x)))
illb <- apply(p3b, 1, function(x)ifelse(y==1, log(x), log(1-x)))

## ------------------------------------------------------------------------
mean(2*(colSums(illa) - colSums(illb)) <  qchisq(0.95, 4))

## ------------------------------------------------------------------------
illD <- illa > illb
summary(colMeans(illD))
mean(colMeans(illD) < .46)

## ------------------------------------------------------------------------
hb3 <- apply(b3b, 2, coda:::HPDinterval.mcmc)
coef3 <- apply(b3b, 2, mean)
p1 <- apply(b3b, 2, function(x)mean(x > 0))
phb3 <- paste("(", apply(round(hb3, 2), 2, paste, collapse=", "), ")", sep="")
pb3b <- round(coef3, 2)
names(pb3b) <- c("Intercept", "Blackmun Distance to Median", "Blackmun's Clerk's Recommendation", "Memo Clerk Recommendation", "Blackmun Distance to Memo Clerk", "Cert-worthiness", "Lower Court Direction", "Memo Clerk Recommendation x Blackmun Distance to Memo Clerk")
shuffle(pb3b, phb3)
coefs[[3]] <- shuffle(pb3b, phb3)

## ------------------------------------------------------------------------
b3bm <- colMeans(b3b)
m3a <- relogit(blackmunvoteoncert2 ~ blackmundistance + blackrecommendation + mcrecommendation*dist_black_mc +
        certworthy1 + lctdir, data=moddat, tau=mean(tau$bv))
p3m <- plogis(model.matrix(m3a) %*% b3bm)
res <- NULL
for(x in seq(0.01, 0.99, length=1000)){
    res <- c(res, mean((p3m > x) ==  model.response(model.frame(m3a))))
}
y3 <- model.response(model.frame(m3a))
my3 <- mean(y3)
pm3 <- max(my3, 1-my3)
pre3 <- (max(res) - pm3)/(1-pm3)

## ------------------------------------------------------------------------
X3ma <- X3mb <- model.matrix(m3a)
X3ma[,6] <- 0
X3mb[,6] <- 1

mes <- plogis(X3mb %*% t(b3b)) - plogis(X3ma %*% t(b3b))
bds <- coda:::HPDinterval.mcmc(colMeans(mes))
bds <- paste("(", round(bds[1], 3), ", ", round(bds[2], 3), ")", sep="")

## ------------------------------------------------------------------------
X3ma <- X3mb <- model.matrix(m3a)
X3ma[,3] <- 0
X3mb[,3] <- 1

mes <- plogis(X3mb %*% t(b3b)) - plogis(X3ma %*% t(b3b))
bds <- coda:::HPDinterval.mcmc(colMeans(mes))
bds <- paste("(", round(bds[1], 3), ", ", round(bds[2], 3), ")", sep="")

## ------------------------------------------------------------------------
cert_probs <- array(dim=c(nrow(cert_data), 1500))
cert_probs[which(cc), ] <- t(bart.post)
cert_probs <- as.data.frame(cert_probs)
cert_probs$docket <- cert_data$docket

## ------------------------------------------------------------------------
cert_post <-  left_join(data.frame(docket=moddat$docket), cert_probs)
cert_post <- cert_post[,-1]
moddat$certworthy1 <- rowMeans(cert_post)

## ------------------------------------------------------------------------
w <- which(complete.cases(moddat) & complete.cases(cert_post))
cert_post <- cert_post[w,] 
moddat <- moddat[w,] 
b1 <- array(dim=c(1500, 4))
set.seed(12294292)
for(i in 1:1500){
    m1a <- relogit(mcrecommendation ~ lctdir + cert_post[,i] + mcdist, data=moddat, tau=tau$mc[i])
    b1[i, ] <- mvrnorm(1, coef(m1a), vcov(m1a))
}

h1 <- apply(b1, 2, coda:::HPDinterval.mcmc)
coef1 <- apply(b1, 2, mean)
p1 <- apply(b1, 2, function(x)mean(x > 0))
ph1 <- paste("(", apply(round(h1, 2), 2, paste, collapse=", "), ")", sep="")
pb1 <- round(coef1, 2)
names(pb1) <- c("Intercept", "Lower Court Direction", "Cert Worthiness", "Distance to Court Median")
coefs[[1]] <- cbind(coefs[[1]], shuffle(pb1, ph1))

## ------------------------------------------------------------------------
b2 <- array(dim=c(1500, 7))
set.seed(692042)
    for(i in 1:1500){
    m2a <- relogit(blackrecommendation ~ blackmundistance + mcrecommendation*dist_black_mc + cert_post[,i] + lctdir, data=moddat, tau=tau$bc[i])
    b2[i, ] <- mvrnorm(1, coef(m2a), vcov(m2a))
}

h2 <- apply(b2, 2, coda:::HPDinterval.mcmc)
coef2 <- apply(b2, 2, mean)
p1 <- apply(b2, 2, function(x)mean(x > 0))
ph2 <- paste("(", apply(round(h2, 2), 2, paste, collapse=", "), ")", sep="")
pb2 <- round(coef2, 2)
names(pb2) <- c("Intercept", "Blackmun Distance to Median", "Memo Clerk Recommendation", "Blackmun Distance to Memo Clerk", "Cert-worthiness", "Lower Court Direction", "Memo Clerk Recommendation x Blackmun Distance to Memo Clerk")
coefs[[2]] <- cbind(coefs[[2]], shuffle(pb2, ph2))

## ------------------------------------------------------------------------
b3b <- array(dim=c(1500, ncol(model.matrix(m3b))))
set.seed(6603)
for(i in 1:1500){
   m3b <- relogit(blackmunvoteoncert2 ~ blackmundistance + blackrecommendation + mcrecommendation*dist_black_mc +
        cert_post[,i] + lctdir, data=moddat, tau=tau$bv[i])
   b3b[i,] <- mvrnorm(1, coef(m3b), vcov(m3b))
}

hb3 <- apply(b3b, 2, coda:::HPDinterval.mcmc)
coef3 <- apply(b3b, 2, mean)
p1 <- apply(b3b, 2, function(x)mean(x > 0))
phb3 <- paste("(", apply(round(hb3, 2), 2, paste, collapse=", "), ")", sep="")
pb3b <- round(coef3, 2)
names(pb3b) <- c("Intercept", "Blackmun Distance to Median", "Blackmun's Clerk's Recommendation", "Memo Clerk Recommendation", "Blackmun Distance to Memo Clerk", "Cert-worthiness", "Lower Court Direction", "Memo Clerk Recommendation x Blackmun Distance to Memo Clerk")
coefs[[3]] <- cbind(coefs[[3]], shuffle(pb3b, phb3))

## ------------------------------------------------------------------------
coefs[[1]]

coefs[[2]]

coefs[[3]]

